// here we have an exogenous rule for export tax

var w1,w2,T1,T2,pi12,pi21,x1,x2,r1,r2,pi11,pi22,P2,riskfree1,riskfree2,taux;

varexo sa1,sa2;

parameters betad, theta1, dbar, sigmaH,sigmaF, alpha2, L1, L2,delta1, rhoa, rho,alpha1,gL,std_ea,rbest,extax_ss,lrss,Ts1,Ts2;

load setpara_ramsey.mat

ncountries     =  mpara(1); % number of countries
nsectors       =  mpara(2);
betad   = mpara(3);     % discount factor
theta1  = mpara(4);
sigmaH = mpara(5);
gL      = mpara(6);
dbar   = mpara(7);
delta1 = mpara(8);
L1     = mpara(9);
L2      = mpara(10);
alpha1 = mpara(11);
alpha2  = mpara(12);
rhoa    = mpara(13);
std_ea  = mpara(14);
extax_ss    = mpara(15);
lrss   = mpara(16);
Ts1    = mpara(17);
Ts2   = mpara(18);
sigmaF = mpara(29);

rho   = betad*(1-delta1);
rbest = delta1/(theta1*(1-betad*(1-delta1))+delta1);



model;


[name = 'P1 definition(1)']
(exp(T1)*exp(w1)^(-theta1)+exp(T2)*(exp(w2)*dbar)^(-theta1))^(-1/theta1)=1;

[name = 'P2 definition(2)']
P2 = -1/theta1*log(exp(T1)*(exp(w1)*(1+taux)*dbar)^(-theta1)+exp(T2)*(exp(w2))^(-theta1));

[name = 'balance of payment(3)']
pi12*exp(x1)=pi21*exp(x2);

[name = 'optimal innovation(4)']
exp(w2)/(exp(sa2)*alpha2) = 1/theta1*exp(w2)*(L2-exp(r2)*L2)/exp(T2)+betad*(1-delta1)*(exp(x2(+1))^(-sigmaF)*exp(P2(+1))^(sigmaF-1)*exp(w2(+1))/(exp(sa2(+1))*alpha2)*exp(x2)^sigmaF*exp(P2)^(1-sigmaF));

[name = 'expenditure x2 (5)']
exp(x2) = (1+theta1)/theta1*exp(w2)*(1-exp(r2))*L2;

[name = 'expenditure x1 (6)']
exp(x1) = (1+theta1)/theta1*exp(w1)*(1-exp(r1))*L1+taux/(1+taux)*pi21*exp(x2);

[name = 'evolution T1(7)']
exp(T1) = exp(T1(-1))*(1-delta1)+exp(sa1)*alpha1*exp(r1)*L1;

[name = 'expenditure T2 (8)']
exp(T2) = exp(T2(-1))*(1-delta1)+exp(sa2)*alpha2*exp(r2)*L2;


//pi11 = exp(T1(-1))*(exp(w1))^(-theta1)/(exp(T1(-1))*exp(w1)^(-theta1)+exp(T2(-1))*(exp(w2)*dbar)^(-theta1));
//pi22 = exp(T2(-1))*(exp(w2))^(-theta1)/(exp(T1(-1))*(exp(w1)*(1+taux)*dbar)^(-theta1)+exp(T2(-1))*(exp(w2))^(-theta1));
//(9, 10, 11, 12)
[name = 'expenditure share (9)-(12)']
pi11 = exp(T1)*(exp(w1))^(-theta1)/(exp(T1)*exp(w1)^(-theta1)+exp(T2)*(exp(w2)*dbar)^(-theta1));
pi22 = exp(T2)*(exp(w2))^(-theta1)/(exp(T1)*(exp(w1)*(1+taux)*dbar)^(-theta1)+exp(T2)*(exp(w2))^(-theta1));
pi12 = 1-pi11;
pi21 = 1-pi22;

[name = 'optimal innovation choice (13)-(14)']
exp(w1)/(alpha1)=1/theta1*exp(w1)*(1-exp(r1))*L1/exp(T1)+betad*(1-delta1)*(exp(x1(+1))^(-sigmaH)*exp(w1(+1))/(alpha1)*exp(x1)^sigmaH);


exp(riskfree1)  = exp(x1)^(-sigmaH)/(betad*exp(x1(+1))^(-sigmaH));
exp(riskfree2)  = exp(x2)^(-sigmaF)*exp(P2)^(1-sigmaF)/(betad*exp(x2(+1))^(-sigmaF)*exp(P2(+1))^(1-sigmaF));
  
//sa1 = rhoa*sa1(-1)+std_ea*ea1;
//sa2 = rhoa*sa2(-1)+std_ea*ea2;

taux = 1/(theta1*pi22);

end;


initval;


T1 = Ts1;
T2 = Ts2;

end;

endval;
sa1 = 0;
sa2 = 0;


r1 = lrss;

end;
steady;


perfect_foresight_setup(periods=2000);
perfect_foresight_solver(stack_solve_algo=6);




